#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# This file is to replicate analyses reported article titled:
# "Does the US Really Embolden Its Allies? Evidence from a Survey Experiment in Japan"

# by Yasuki Kudo and Nguyen Cao Viet Hung.

# Last update: 2025-07-05 by Yasuki Kudo
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

rm(list=ls())

# Packages
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

#devtools::install_github("JaehyunSong/BalanceR", force = TRUE)
library(BalanceR)
library(ggrepel) 
library(ggpp)
library(ggpubr)
library(marginaleffects)
library(modelsummary)
library(tidyverse)

# Data
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

load("data/kudo_nguyen_caoJCR2025_data.RData")


# Analysis 1
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

# Visualization (Figure 2)

plot.policy.approval <- df |>
  filter(attention == "pass") |>
  group_by(t1) |>
  summarise(p = mean(policy_approve), n = n()) |>
  ungroup() |>
  mutate(se = sqrt((p * (1 - p)) / n)) |>
  ggplot(aes(
    x     = t1,
    y     = p * 100,
    ymax  = (p + se * 1.96) * 100,
    ymin  = (p - se * 1.96) * 100,
    label = round(p * 100, 1)
  )) +
  theme_bw() +
  geom_pointrange(linewidth = 1, size = 0.8) +
  geom_text(aes(hjust = -0.5), size = 4) +
  scale_x_discrete(labels = c(
    "control"      = "C",
    "verbal"       = "C + \n Verbal",
    "verbal_force" = "C + \n Verbal + Force"
  )) +
  theme(
    text          = element_text(size = 10, color = "black"),
    axis.text.x   = element_text(size = 10, color = "black"),
    axis.text.y   = element_text(size = 8, color = "black"),
    plot.title    = element_text(size = 10)
  ) +
  ylim(0, 100) +
  ylab("Approval for Dispatching SDF") +
  xlab("Treatment Groups in First Stage \n (US Signals)")

print(plot.policy.approval)

ggsave("figure/figure2.pdf", plot = plot.policy.approval, width  = 6,
       height = 5, units = "in", dpi = 2000)

ggsave("figure/figure2.png", plot = plot.policy.approval, width  = 6,
       height = 5, units = "in", dpi = 2000)

# ANOVA (footnote 22)
#---------------------------------------------------------------------------------------#

aov(policy_support ~ t1, data = df, subset = attention == "pass") |> summary()

# Regression analyses (Appendix 6.1)
#---------------------------------------------------------------------------------------#

# control vars
control <- c(
  "as.factor(male)", "as.factor(cdegree)", "as.factor(income)",
  "as.factor(region)", "know_n", "news_pol", "as.factor(ldp_support)", "age"
  )

### Logit
# no controls
reg.pol.app.all.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass"
  )

reg.pol.app.all.1.attention <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df
  )

reg.pol.app.hawk.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & hawk_dove == "hawk"
  )

reg.pol.app.neu.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & hawk_dove == "neutral"
  )

reg.pol.app.dove.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & hawk_dove == "dove"
  )

reg.pol.app.sq.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & exp_sdf %in% c(3, 4)
  )

reg.pol.app.agg.1 <- glm(
  policy_approve ~ t1, family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

# controls included
reg.pol.app.all.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass"
  )

reg.pol.app.all.2.attention <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df
  )

reg.pol.app.hawk.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass" & hawk_dove == "hawk"
  )

reg.pol.app.neu.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass" & hawk_dove == "neutral"
  )

reg.pol.app.dove.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass" & hawk_dove == "dove"
  )

reg.pol.app.sq.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass" & exp_sdf %in% c(3, 4)
  )

reg.pol.app.agg.2 <- glm(
  update(policy_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"), data = df,
  subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

### OLS

reg.pol.app.all.1.ols <- lm(
  policy_support ~ t1, data = df,
  subset = attention == "pass"
  )

reg.pol.app.all.2.ols <- lm(
  update(policy_support ~ t1, reformulate(c(".", control))),
  data = df, subset = attention == "pass"
  )

### Tables

var.name <- c(
  t1verbal                         = "C + Verbal",
  t1verbal_force                   = "C + Verbal+Force",
  t2develop                        = "Economic Development Rhetoric",
  t2peace                          = "Peaceful People Rhetoric",
  t2sanction                       = "Economic Sanction Rhetoric",
  t2un                             = "UN Mediation Rhetoric",
  `t1verbal:t2develop`             = "C + Verbal x Economic Development ",
  `t1verbal_force:t2develop`       = "C + Verbal+Force x Economic Development ",
  `t1verbal:t2peace`               = "C + Verbal x Peaceful People",
  `t1verbal_force:t2peace`         = "C + Verbal+Force x Peaceful People",
  `t1verbal:t2sanction`            = "C + Verbal x Economic Sanction",
  `t1verbal_force:t2sanction`      = "C + Verbal+Force x Economic Sanction",
  `t1verbal:t2un`                  = "C + Verbal x UN Mediation",
  `t1verbal_force:t2un`            = "C + Verbal+Force x UN Mediation",
  `as.factor(male)`                = "Gender (male or not)",
  `as.factor(cdegree)`             = "College Degree",
  `as.factor(income)2`             = "Income Level 2",
  `as.factor(income)3`             = "Income Level 3",
  `as.factor(income)4`             = "Income Level 4",
  ldp_support                      = "LDP Support",
  `as.factor(region)tohoku`        = "Region = Tohoku",
  `as.factor(region)kanto`         = "Region = Kanto",
  `as.factor(region)chubu`         = "Region = Chubu",
  `as.factor(region)kinki`         = "Region = Kinki",
  `as.factor(region)chugoku_shikoku` = "Region = Chugoku & Shikoku",
  `as.factor(region)kyushu_okinawa`  = "Region = Kyushu & Okinawa",
  know_n                           = "Knowledge",
  news_pol                         = "Interest in politics",
  `as.factor(ldp_support)1`        = "PID",
  age                              = "Age"
  )

# Table A1
modelsummary(
  models = list(
    `Logit`            = reg.pol.app.all.1,
    `Logit`            = reg.pol.app.all.2,
    `Logit (All res.)` = reg.pol.app.all.1.attention,
    `Logit (All res.)` = reg.pol.app.all.2.attention,
    `OLS`              = reg.pol.app.all.1.ols,
    `OLS`              = reg.pol.app.all.2.ols
  ),
  gof_map = c("nobs"), coef_map = var.name, stars = c('*' = .05), 
  align = 'ldddddd', title = "Table A1"
  )

# Table A2
modelsummary(
  models = list(
    `Hawk`             = reg.pol.app.hawk.1,
    `Hawk`             = reg.pol.app.hawk.2,
    `Neutral`          = reg.pol.app.neu.1,
    `Neutral`          = reg.pol.app.neu.2,
    `Dove`             = reg.pol.app.dove.1,
    `Dove`             = reg.pol.app.dove.2,
    `Less Escalatory`  = reg.pol.app.sq.1,
    `Less Escalatory`  = reg.pol.app.sq.2,
    `More Escalatory`  = reg.pol.app.agg.1,
    `More Escalatory`  = reg.pol.app.agg.2
  ),
  gof_map  = c("nobs"), coef_map = var.name, stars = c('*' = .05),
  align = 'ldddddddddd', title = "Table A2"
  )

# Analysis 2 (1): US commitments and approval for back down 
# (Figure 2 and related analyses)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

# ANOVA (footnote 27)
#---------------------------------------------------------------------------------------#

aov(bd_support ~ t1, data = df, subset = attention == "pass" & t2 == "control") |> summary()

# Regression Analyses 
#---------------------------------------------------------------------------------------#

### Logit

mod.bd.nc <- glm(
  bd_approve ~ t1,
  family = binomial(link = "logit"),
  data   = df, subset = (t2 == "control" & attention == "pass")
  )

mod.bd.c <- glm(
  update(bd_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data   = df, subset = (t2 == "control" & attention == "pass")
  )

mod.bd.nc.attention <- glm(
  bd_approve ~ t1,
  family = binomial(link = "logit"),
  data   = df, subset = t2 == "control"
  )

mod.bd.c.attention <- glm(
  update(bd_approve ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data   = df, subset = t2 == "control"
  )

### OLS

mod.bd.nc.ols <- lm(
  bd_support ~ t1,
  data   = df, subset = (t2 == "control" & attention == "pass")
  )

mod.bd.c.ols <- lm(
  update(bd_support ~ t1, reformulate(c(".", control))),
  data   = df, subset = (t2 == "control" & attention == "pass")
  )

### Figure 3

comp.bd.nc <- comparisons(mod.bd.nc, variable = "t1", newdata = datagrid()) |>
  dplyr::select(contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "not included")

comp.bd.c  <- comparisons(mod.bd.c,  variable = "t1", newdata = datagrid()) |>
  dplyr::select(contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "included")

rbind(comp.bd.nc, comp.bd.c) |>
  mutate(
    x = ifelse(contrast == "verbal - control", "C + Verbal", "C + Verbal + Force")
  ) |>
  ggplot(aes(x, estimate * 100, shape = Control, fill = Control)) +
  theme_bw() +
  geom_pointrange(aes(ymin = conf.low * 100, ymax = conf.high * 100),
                  position = position_dodge(0.5), size = 1, linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  scale_x_discrete(
    limits = c("verbal_force" = "C + Verbal + Force", "verbal" = "C + Verbal")
  ) +
  scale_shape_manual(values = c(16, 21)) +
  scale_fill_manual(values = c("black", "white")) +
  theme(
    text            = element_text(size = 10),
    axis.text.y     = element_text(size = 10, color = "black"),
    axis.text.x     = element_text(size = 8, color = "black"),
    legend.position = "top"
  ) +
  coord_flip() +
  guides(shape = guide_legend("Control variables"), fill = guide_legend("Control variables")) +
  ylab("Average Effects on Percentage Approval (% point)\n(Baseline = C)") +
  xlab("Treatment Groups in First Stage\n(US Signals)") ->
  bd.effect.plot

print(bd.effect.plot)

ggsave("figure/figure4.pdf", bd.effect.plot, 
       width = 6, height = 4, units = "in", dpi = 2000)

ggsave("figure/figure4.png", bd.effect.plot, 
       width = 6, height = 4, units = "in", dpi = 2000)

### Table A3

modelsummary(
  models  = list(
    `Logit`            = mod.bd.nc,
    `Logit`            = mod.bd.c,
    `Logit (All res.)` = mod.bd.nc.attention,
    `Logit (All res.)` = mod.bd.c.attention,
    `OLS`              = mod.bd.nc.ols,
    `OLS`              = mod.bd.c.ols
  ),
  gof_map = "nobs", coef_map = var.name,
  stars = c('*' = .05), align  = 'ldddddd', title = "Table A3"
  )

# Analysis 2 (2): US commitments and side-stepping strategy
# (Figure 4 and related analyses)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

# Regression analyses
#---------------------------------------------------------------------------------------#

### Logit

mod.side.nc <- glm(
  bd_approve ~ t1 * t2,
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass"
  )

mod.side.c <- glm(
  update(bd_approve ~ t1 * t2, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass"
  )

mod.side.nc.attention <- glm(
  bd_approve ~ t1 * t2,
  family = binomial(link = "logit"),
  data = df
  )

mod.side.c.attention <- glm(
  update(bd_approve ~ t1 * t2, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data = df
  )

### OLS

mod.side.nc.ols <- lm(
  bd_support ~ t1 * t2,
  data = df, subset = attention == "pass"
  )

mod.side.c.ols <- lm(
  update(bd_support ~ t1 * t2, reformulate(c(".", control))),
  data = df, subset = attention == "pass"
  )

### Figure 4

comp.side.nc <- comparisons(
  model     = mod.side.nc,
  variables = "t2",
  newdata   = datagrid(t1 = c("control", "verbal", "verbal_force", t2 = "control"))
  ) |>
  dplyr::select(
    "t1", "contrast", "estimate", "std.error",
    "p.value", "conf.low", "conf.high"
  ) |>
  mutate(Control = "not included")

comp.side.c <- comparisons(
  model     = mod.side.c,
  variables = "t2",
  newdata   = datagrid(t1 = c("control", "verbal", "verbal_force", t2 = "control"))
  ) |>
  dplyr::select(
    "t1", "contrast", "estimate", "std.error",
    "p.value", "conf.low", "conf.high"
  ) |>
  mutate(Control = "included")

rbind(comp.side.nc, comp.side.c) |>
  mutate(
    contrast_nme = case_when(
      contrast == "develop - control"  ~ "Economic Development",
      contrast == "peace - control"    ~ "Peaceful People",
      contrast == "sanction - control" ~ "Economic Sanction",
      TRUE                             ~ "UN Mediation"
    )
  ) |>
  ggplot(aes(
    x     = t1,
    y     = estimate * 100,
    shape = Control,
    fill  = Control
  )) +
  theme_bw() +
  geom_pointrange(
    aes(
      ymax = conf.high * 100,
      ymin = conf.low  * 100
    ),
    position  = position_dodge(width = 0.5),
    size      = 1,
    linewidth = 1
  ) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  scale_x_discrete(
    labels = c(
      "control"      = "C",
      "verbal"       = "C + Verbal",
      "verbal_force" = "C + Verbal + Force"
    ),
    limits = c("verbal_force", "verbal", "control")
  ) +
  scale_shape_manual(values = c(16, 21)) +
  scale_fill_manual(values  = c("black", "white")) +
  theme(
    text            = element_text(size = 10, color = "black"),
    axis.text.y     = element_text(size = 10, color = "black"),
    axis.text.x     = element_text(size = 8,  color = "black"),
    plot.title      = element_text(size = 10),
    legend.position = "top"
  ) +
  facet_wrap(~contrast_nme) +
  guides(
    shape = guide_legend(title = "Control variables"),
    fill  = guide_legend(title = "Control variables")
  ) +
  coord_flip() +
  ylab(
    "Average Effects on Percentage Approval (% point) \n (Baseline = Backdown w/o Any Sidestepping Strategies)"
  ) +
  xlab(
    "Treatment Groups in First Stage \n (US Signals)"
  ) ->
  plot.sidestep.support

print(plot.sidestep.support)

ggsave("figure/figure5.pdf", plot = plot.sidestep.support,
  width = 6, height = 6, units = "in", dpi = 2000)

ggsave("figure/figure5.png", plot = plot.sidestep.support,
       width = 6, height = 6, units = "in", dpi = 2000)

### Table A4

modelsummary(
  models = list(
    `Logit`            = mod.side.nc,
    `Logit`            = mod.side.c,
    `Logit (All res.)` = mod.side.nc.attention,
    `Logit (All res.)` = mod.side.c.attention,
    `OLS`              = mod.side.nc.ols,
    `OLS`              = mod.side.c.ols
  ),
  gof_map  = c("nobs"), coef_map = var.name,
  stars = c('*' = .05), align = "ldddddd", title = "Table A4"
  )

# Additional Analyses
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

# Association between the evaluation of power balance and the treatments 
# (footnote 24) 
#---------------------------------------------------------------------------------------#

aov(powbal ~ t1, data = df, subset = attention == "pass") |> summary()
aov(powbal ~ t2, data = df, subset = attention == "pass") |> summary()


# US military intervention (Figure 3 and related analyses) 
#---------------------------------------------------------------------------------------#

### Average perceived likelihood of US military intervention

subset(df, attention == "pass" & t1 == "control")$lik_milint_bi |> summary() # binary
subset(df, attention == "pass" & t1 == "control")$lik_milint |> summary()   # continuous

### Regression analyses

# Logit

mod.milint.nc <- glm(
  lik_milint_bi ~ t1,
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass"
  )

mod.milint.c <- glm(
  update(lik_milint_bi ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass"
  )

mod.milint.aggressive.nc <- glm(
  lik_milint_bi ~ t1,
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

mod.milint.aggressive.c <- glm(
  update(lik_milint_bi ~ t1, reformulate(c(".", control))),
  family = binomial(link = "logit"),
  data = df, subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

# OLS

mod.milint.nc.ols <- lm(
  lik_milint ~ t1,
  data = df, subset = attention == "pass"
  )

mod.milint.c.ols <- lm(
  update(lik_milint ~ t1, reformulate(c(".", control))),
  data = df, subset = attention == "pass"
)

mod.milint.aggressive.nc.ols <- lm(
  lik_milint ~ t1,
  data = df, subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

mod.milint.aggressive.c.ols <- lm(
  update(lik_milint ~ t1, reformulate(c(".", control))),
  data = df, subset = attention == "pass" & exp_sdf %in% c(1, 2)
  )

### Figure 2

comp.milint.nc <- comparisons(mod.milint.nc, variable = "t1", newdata = datagrid()) |>
  dplyr::select(t1, contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "not included")

comp.milint.c  <- comparisons(mod.milint.c,  variable = "t1", newdata = datagrid()) |>
  dplyr::select(t1, contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "included")

rbind(comp.milint.nc, comp.milint.c) |>
  mutate(x = ifelse(contrast == "verbal - control", "C + Verbal", "C + Verbal + Force")) |>
  ggplot(aes(x, estimate * 100, shape = Control, fill = Control)) +
    geom_pointrange(aes(ymin = conf.low * 100, ymax = conf.high * 100),
                    position = position_dodge(0.5),
                    size = 1,
                    linewidth = 1) +
    geom_hline(yintercept = 0, linetype = "dotted") +
    scale_shape_manual(values = c(16, 21)) +
    scale_fill_manual(values = c("black", "white")) +
    scale_x_discrete(limits = c("C + Verbal + Force", "C + Verbal")) +
    coord_flip() +
    theme_bw() +
    theme(
      text        = element_text(size = 10),
      axis.text.x = element_text(size = 8, color = "black"),
      axis.text.y = element_text(size = 10, color = "black"),
      legend.position = "top"
    ) +
    guides(shape = guide_legend("Control variables"), fill = guide_legend("Control variables")) +
    labs(
      x   = "Treatment Groups in First Stage\n(US Signals)",
      y   = "Average Effects on Percentage of Respondents\nwith Score > 50 (% point)"
    ) ->
  plot.us.support

print(plot.us.support)

ggsave("figure/figure3.pdf", plot = plot.us.support,
       width = 6, height = 4, units = "in", dpi = 2000)

ggsave("figure/figure3.png", plot = plot.us.support,
       width = 6, height = 4, units = "in", dpi = 2000)

### Figure A6

comp.milint.aggressive.nc <- comparisons(mod.milint.aggressive.nc, variable = "t1", newdata = datagrid()) |>
  dplyr::select(t1, contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "not included")

comp.milint.aggressive.c <- comparisons(mod.milint.aggressive.c, variable = "t1", newdata = datagrid()) |>
  dplyr::select(t1, contrast, estimate, std.error, p.value, conf.low, conf.high) |>
  mutate(Control = "included")

rbind(comp.milint.aggressive.nc, comp.milint.aggressive.c) |>
  mutate(x = ifelse(contrast == "verbal - control", "C + Verbal", "C + Verbal + Force")) |>
  ggplot(aes(x, estimate * 100, shape = Control, fill = Control)) +
    geom_pointrange(aes(ymin = conf.low * 100, ymax = conf.high * 100),
                    position = position_dodge(0.5), 
                    size = 1,
                    linewidth = 1) +
    geom_hline(yintercept = 0, linetype = "dotted") +
    scale_shape_manual(values = c(16, 21)) +
    scale_fill_manual(values  = c("black", "white")) +
    scale_x_discrete(limits   = c("C + Verbal + Force", "C + Verbal")) +
    coord_flip() +
    theme_bw() +
    theme(
      text        = element_text(size = 10),
      axis.text.x = element_text(size = 8, color = "black"),
      axis.text.y = element_text(size = 10, color = "black"),
      legend.position = "top"
    ) +
    guides(shape = guide_legend("Control variables"), fill = guide_legend("Control variables")) +
    labs(
      x   = "Treatment Groups in First Stage\n(US Signals)",
      y   = "Average Effects on Percentage of Respondents\nwith Score > 50 (% point)"
    ) ->
  plot.us.support.aggressive

print(plot.us.support.aggressive)

ggsave("figure/figureA6.pdf", plot = plot.us.support.aggressive, 
       width = 6, height = 4, units = "in", dpi = 2000)

ggsave("figure/figureA6.png", plot = plot.us.support.aggressive, 
       width = 6, height = 4, units = "in", dpi = 2000)

### Table A5

modelsummary(
  models = list(
    `Logit` = mod.milint.nc,
    `Logit` = mod.milint.c,
    `OLS`   = mod.milint.nc.ols,
    `OLS`   = mod.milint.c.ols,
    `Logit` = mod.milint.aggressive.nc,
    `Logit` = mod.milint.aggressive.c,
    `OLS`   = mod.milint.aggressive.nc.ols,
    `OLS`   = mod.milint.aggressive.c.ols
  ),
  gof_map = "nobs", coef_map = var.name, stars = c('*' = .05), 
  align = 'ldddddddd', title = "Table A5"
  )

# ANOVA for likelihood of a militarized dispute (footnote 25)
#---------------------------------------------------------------------------------------#

aov(lik_mildisp ~ t1, data = df, subset = attention == "pass") |> summary()

# Expected level of SDF's activities (Appendix A3)
#---------------------------------------------------------------------------------------#

### Figure A3

plot.sdf.action <- df |>
  filter(attention == "pass") |>
  ggplot(aes(x = reorder(as.factor(exp_sdf), desc(exp_sdf)))) +
  theme_bw() +
  geom_bar() +
  scale_x_discrete(labels = c(
    "1"  = "Board and seize",
    "2"  = "Obstruction \n using weapons",
    "3"  = "Obstruction \n w/o using weapons",
    "4"  = "Surveillance",
    "5"  = "Others",
    "99" = "Don't know"
  )) +
  theme(
    text           = element_text(size = 10, color = "black"),
    axis.text.x    = element_text(size = 8, color = "black"),
    axis.text.y    = element_text(size = 10, color = "black"),
    plot.title     = element_text(size = 10)
  ) +
  xlab("") +
  ylab("Frequency") +
  coord_flip()

print(plot.sdf.action)

ggsave(
  "figure/figureA3.pdf", plot = plot.sdf.action,
  width = 6, height = 4, units = "in", dpi = 2000
)

ggsave(
  "figure/figureA3.png", plot = plot.sdf.action,
  width = 6, height = 4, units = "in", dpi = 2000
)

### Association between T1/T2 and respondents' expectation of SDF's actions 
# (footnote 5 in Appendix)
df_pass_sdf <- subset(df, attention == "pass" & exp_sdf %in% 1:4)
with(df_pass_sdf, chisq.test(t1, as.factor(exp_sdf)))
with(df_pass_sdf, chisq.test(t2, as.factor(exp_sdf)))


# Balance checks (Appendix A3)
#---------------------------------------------------------------------------------------#
### Figure A2

# Treatment Set 1 (US Signals)
balancet1 <- df |>
  filter(attention == "pass") |>
  mutate(t1_name = case_when(
    t1 == "control" ~ "Control",
    t1 == "verbal"  ~ "Verbal",
    TRUE            ~ "Verbal+Force"
  )) |>
  BalanceR(
    group = t1_name,
    cov   = c(
      Gender          = male,
      Age             = age,
      College_Degree  = cdegree,
      Income_1        = income1,
      Income_2        = income2,
      Income_3        = income3,
      Income_4        = income4,
      Income_5        = income5,
      PID             = ldp_support,
      Hokkaido        = hokkaido,
      Tohoku          = tohoku,
      Kanto           = kanto,
      Chubu           = chubu,
      Kinki           = kinki,
      Chugoku_Shikoku = chugoku_shikoku,
      Kyushu_Okinawa  = kyushu_okinawa
    )
  ) |>
  plot(point.size = 2, text.size = 8) +
  theme(
    legend.position = "bottom",
    text            = element_text(size = 8),
    plot.title      = element_text(size = 8)
  ) +
  guides(col = guide_legend(ncol = 3)) +
  ggtitle("Balance Check for Treatment Set 1 (US Signals)")


# Treatment Set 2 (Side‐stepping strategies)
balancet2 <- df |>
  filter(attention == "pass") |>
  mutate(t2_name = case_when(
    t2 == "control"  ~ "Control",
    t2 == "develop"  ~ "Econ. Devel.",
    t2 == "peace"    ~ "Peaceful People",
    t2 == "sanction" ~ "Econ. Sanc.",
    TRUE             ~ "UN Mediation"
  )) |>
  BalanceR(
    group = t2_name,
    cov   = c(
      Gender          = male,
      Age             = age,
      College_Degree  = cdegree,
      Income_1        = income1,
      Income_2        = income2,
      Income_3        = income3,
      Income_4        = income4,
      Income_5        = income5,
      PID             = ldp_support,
      Hokkaido        = hokkaido,
      Tohoku          = tohoku,
      Kanto           = kanto,
      Chubu           = chubu,
      Kinki           = kinki,
      Chugoku_Shikoku = chugoku_shikoku,
      Kyushu_Okinawa  = kyushu_okinawa
    )
  ) |>
  plot(point.size = 2, text.size = 8) +
  theme(
    legend.position = "bottom",
    text            = element_text(size = 8),
    plot.title      = element_text(size = 8)
  ) +
  guides(col = guide_legend(ncol = 3)) +
  ggtitle("Balance Check for Treatment Set 2 (Sidestepping Strategies)")


balance_plot <- ggarrange(balancet1, balancet2, ncol = 1)

print(balance_plot)

ggsave("figure/figureA2.pdf", plot = balance_plot,
  width = 6, height = 8, units = "in", dpi = 2000)

ggsave("figure/figureA2.png", plot = balance_plot,
       width = 6, height = 8, units = "in", dpi = 2000)


# Sample information (Appendix A6.1)
#---------------------------------------------------------------------------------------#

# Figure A5

plot.hawkdove <- df |>
  ggplot(aes(x = hawk_dove_num)) +
  theme_bw() +
  geom_bar() +
  theme(
    text         = element_text(size = 10, color = "black"),
    axis.text.x  = element_text(size = 8, color = "black"),
    axis.text.y  = element_text(size = 10, color = "black"),
    plot.title   = element_text(size = 10)
  ) +
  xlab("Hawk (1) - Dove (5) Scale") +
  ylab("Frequency")

print(plot.hawkdove)

ggsave("figure/figureA5.pdf", plot = plot.hawkdove,
  width = 5, height = 3.5, units = "in", dpi = 2000)

ggsave("figure/figureA5.png", plot = plot.hawkdove,
       width = 5, height = 3.5, units = "in", dpi = 2000)

# TOST (Appendix A8)
#---------------------------------------------------------------------------------------#

### Figure A7

tost.h1.2 <- function(fit){
  
  m <- seq(0, 0.2, 0.005)
  length.m <- length(m)
  
  for(i in 1:length.m){
    
    if(i == 1) comp.tost <- tibble()
    
    m.i <- m[i]  
    sub.comp <- comparisons(model = fit, variables = "t1", newdata = datagrid(t1 = "control")) |> 
      hypotheses(equivalence = c(-m.i, m.i), conf_level = 0.9) |> 
      dplyr::select(term, contrast, estimate, conf.low, conf.high, p.value.equiv)
    
    comp.tost <- rbind(comp.tost, sub.comp)
    
  }
  
  comp.tost$m <- sort(c(m, m))
  comp.tost ->> comp.tost
  
}

tost.h1.nc <- tost.h1.2(fit = reg.pol.app.all.1) |> mutate(control = "not included")
tost.h1.c <- tost.h1.2(fit = reg.pol.app.all.2) |> mutate(control = "included")

tost.h1.plot <- rbind(tost.h1.nc, tost.h1.c) |>
  mutate(
    x = ifelse(
      contrast == "verbal - control",
      "C + Verbal",
      "C + Verbal + Force"
    )
  ) |>
  ggplot(aes(
    x        = m * 100,
    y        = p.value.equiv,
    color    = x,
    linetype = control
  )) +
  theme_bw() +
  geom_hline(yintercept = 0.05, color = "gray", linewidth = 0.7) +
  geom_line(linewidth = 0.7) +
  scale_y_continuous(limits = c(0, 1)) +
  theme(
    text             = element_text(size = 10, color = "black"),
    axis.text        = element_text(size = 8, color = "black"),
    plot.title       = element_text(size = 10),
    legend.position  = "right",
    legend.box       = "vertical",
    legend.spacing.y = unit(0, "mm"),
    panel.grid       = element_blank()
  ) +
  guides(
    color    = guide_legend(title = "US signals"),
    linetype = guide_legend(title = "Control variables")
  ) +
  labs(
    x = "Equivalence value (m)\n(Hypothetical treatment effect of US signals (% point))",
    y = "TOST p-value"
  )

print(tost.h1.plot)

ggsave("figure/figureA7.pdf", plot = tost.h1.plot,
       width = 6, height = 4, units = "in", dpi = 2000)

ggsave("figure/figureA7.png", plot = tost.h1.plot,
       width = 6, height = 4, units = "in", dpi = 2000)  

### Figure A8
tost.h2.nc <- tost.h1.2(fit = mod.bd.nc) |> mutate(control = "not included")
tost.h2.c <- tost.h1.2(fit = mod.bd.c) |> mutate(control = "included")

tost.h2.plot <- rbind(tost.h2.c, tost.h2.nc) |>
  mutate(
    x = ifelse(
      contrast == "verbal - control",
      "C + Verbal",
      "C + Verbal + Force"
    )
  ) |>
  ggplot(aes(
    x        = m * 100,
    y        = p.value.equiv,
    color    = x,
    linetype = control
  )) +
  theme_bw() +
  geom_hline(
    yintercept = 0.05,
    color      = "gray",
    linewidth  = 0.7
  ) +
  geom_line(linewidth = 0.7) +
  theme(
    text             = element_text(size = 10, color = "black"),
    axis.text.y      = element_text(size = 8, color = "black"),
    axis.text.x      = element_text(size = 8, color = "black"),
    plot.title       = element_text(size = 10),
    legend.position  = "right",
    legend.box       = "vertical",
    legend.spacing.y = unit(0, "mm"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  ) +
  guides(
    color    = guide_legend(title = "US signals"),
    linetype = guide_legend(title = "Control variables")
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    x = "Equivalence value (m)\n(Hypothetical treatment effect of US signals (% point))",
    y = "TOST p-value"
  )
print(tost.h2.plot)

ggsave("figure/figureA8.pdf", plot = tost.h2.plot, 
       width = 6, height = 4, units = "in", dpi = 2000)

ggsave("figure/figureA8.png", plot = tost.h2.plot, 
       width = 6, height = 4, units = "in", dpi = 2000)  

### Figure A9
tost.h3 <- function(fit){
  
  sstep.list <- rep(c("develop", "peace", "un", "sanction"), 2)
  signal.list <- c(rep("verbal", 4), rep("verbal_force", 4))
  m <- seq(0, 0.2, 0.005)
  length.m <- length(m)
  
  for(j in 1:8){
    
    if(j == 1) comp.tost.pooled <- tibble()
    
    sstep <- sstep.list[j]
    signal <- signal.list[j]
    
    for(i in 1:length.m){
      
      if(i == 1) comp.tost <- tibble()
      
      m.i <- m[i]  
      
      sub.comp <- comparisons(model = fit, variables = list("t2" = c(sstep, "control")), 
                              newdata = datagrid(t2 = "control", t1 = c(signal, "control")), hypothesis = "pairwise") |> 
        hypotheses(equivalence = c(-m.i, m.i), conf_level = 0.9) |> 
        dplyr::select(term, estimate, conf.low, conf.high, p.value.equiv) |> 
        mutate(sstep = sstep,
               signal = signal)
      
      comp.tost <- rbind(comp.tost, sub.comp)
      
    }
    
    comp.tost$m <- m
    comp.tost.pooled <- rbind(comp.tost, comp.tost.pooled)
    cat(j) 
    
  }
  
  comp.tost.pooled ->> comp.tost.pooled
  
}

tost.h3.c <- tost.h3(fit = mod.side.c) |> mutate(control = "included")
tost.h3.nc <- tost.h3(fit = mod.side.nc)|> mutate(control = "not included")

tost.h3.plot <- rbind(tost.h3.c, tost.h3.nc) |>
  mutate(
    signal = ifelse(
      signal == "verbal",
      "from C to C + Verbal",
      "from C to C + Verbal + Force"
    ),
    sstep = case_when(
      sstep == "control" ~ "Control",
      sstep == "develop" ~ "Economic Development",
      sstep == "peace"   ~ "Peaceful People",
      sstep == "sanction"~ "Economic Sanction",
      TRUE               ~ "UN Mediation"
    )
  ) |>
  ggplot(aes(
    x        = m * 100,
    y        = p.value.equiv,
    linetype = control,
    color    = signal
  )) +
  theme_bw() +
  geom_hline(yintercept = 0.05, color = "gray", linewidth = 0.7) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~sstep) +
  scale_y_continuous(limits = c(0, 1)) +
  theme(
    text               = element_text(size = 10, color = "black"),
    axis.text          = element_text(size = 8, color = "black"),
    plot.title         = element_text(size = 10),
    panel.grid         = element_blank(),
    legend.position    = "top",
    legend.box         = "vertical",
    legend.spacing.y    = unit(-3, "mm")
  ) +
  guides(
    color    = guide_legend(title = "Shifts in US signals"),
    linetype = guide_legend(title = "Control variables")
  ) +
  labs(
    x = "Equivalence value (m)\n(Hypothetical shift in treatment effect of sidestepping strategy\nalong with a shift in US signals (% point))",
    y = "TOST p-value"
  )

print(tost.h3.plot)

ggsave("figure/figureA9.pdf", plot = tost.h3.plot, 
       width = 6, height = 6, units = "in", dpi = 2000)  

ggsave("figure/figureA9.png", plot = tost.h3.plot, 
       width = 6, height = 6, units = "in", dpi = 2000)  


# Follow-up survey (Appendix A5)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

### Proportions

# proportion of "escalate" and "escalate somewhat"
p.esc <- ((table(df_followup$escalate)[4] + table(df_followup$escalate)[5])/nrow(df_followup))
se.esc <- sqrt(p.esc*(1-p.esc)/nrow(df_followup))
paste(round(p.esc-1.96*se.esc, 6), round(p.esc, 6), round(p.esc+1.96*se.esc, 6))

# proportion of "increase" and "increase somewhat"
p.disp <- (table(df_followup$mildisp)[4] + table(df_followup$mildisp)[5])/nrow(df_followup)
se.disp <- sqrt(p.disp*(1-p.disp)/nrow(df_followup))
paste(round(p.disp-1.96*se.disp, 6), round(p.disp, 6), round(p.disp+1.96*se.disp, 6))

# both
p.both <- nrow(df_followup[df_followup$escalate %in% c(4, 5) | df_followup$mildisp %in% c(4, 5),])/nrow(df_followup)
se.both <- sqrt(p.both*(1-p.both)/nrow(df_followup))
paste(round(p.both-1.96*se.both, 6), round(p.both, 6), round(p.both+1.96*se.both, 6))

### Figure A4

df_followup.esc <- df_followup |> 
  (\(x) mutate(x, n = nrow(x)))() |> 
  group_by(escalate, escalate_name) |>
  reframe(p = n(),
          n = mean(n)) |> 
  mutate(x = p/n,
         y = "Will PM's decision escalate or de-escalate \n the conflict?") |> 
  rename(label = escalate_name,
         grp = escalate)

df_followup.disp <- df_followup |> 
  (\(x) mutate(x, n = nrow(x)))() |> 
  group_by(mildisp, mildisp_name) |>
  reframe(p = n(),
          n = mean(n)) |> 
  mutate(x = p/n,
         y = "Will PM's decision increase or decrease \n the likelihood of a military clash?") |> 
  rename(label = mildisp_name,
         grp = mildisp)

plot.followup.bar <- rbind(df_followup.esc, df_followup.disp) |>
  ggplot(aes(x, y, group = as.factor(grp))) +
  theme_bw() +
  geom_col(aes(fill = as.factor(grp)), width = 0.3) +
  geom_label_repel(
    aes(label = label),
    position     = position_stacknudge(vjust = 0.5, y = 0.5),
    size         = 3.7,
    label.padding = 0.2,
    label.size   = NA
  ) +
  scale_fill_discrete(guide = "none") +
  scale_x_continuous(breaks = seq(0, 1, 0.2)) +
  theme(
    text          = element_text(size = 12, color = "black"),
    axis.text.x   = element_text(size = 10, color = "black"),
    axis.text.y   = element_text(size = 12, color = "black"),
    plot.title    = element_text(size = 12)
  ) +
  ylab("") +
  xlab("Proportions")

print(plot.followup.bar)

ggsave("figure/figureA4.pdf", plot = plot.followup.bar, 
       width = 10, height = 4, units = "in", dpi = 2000)

ggsave("figure/figureA4.png", plot = plot.followup.bar, 
       width = 10, height = 4, units = "in", dpi = 2000)

# [End of the script]  
